Modeling time-varying phytoplankton subsidy reveals at-risk species in a Chilean intertidal ecosystem

The allometric trophic network (ATN) framework for modeling population dynamics has provided numerous insights into ecosystem functioning in recent years. Herein we extend ATN modeling of the intertidal ecosystem off central Chile to include empirical data on pelagic chlorophyll-a concentration. This intertidal community requires subsidy of primary productivity to support its rich ecosystem. Previous work models this subsidy using a constant rate of phytoplankton input to the system. However, data shows pelagic subsidies exhibit highly variable, pulse-like behavior. The primary contribution of our work is incorporating this variable input into ATN modeling to simulate how this ecosystem may respond to pulses of pelagic phytoplankton. Our model results show that: (1) closely related sea snails respond differently to phytoplankton variability, which is explained by the underlying network structure of the food web; (2) increasing the rate of pelagic-intertidal mixing increases fluctuations in species’ biomasses that may increase the risk of local extirpation; (3) predators are the most sensitive species to phytoplankton biomass fluctuations, putting these species at greater risk of extirpation than others. Finally, our work provides a straightforward way to incorporate empirical, time-series data into the ATN framework that will expand this powerful methodology to new applications.

of this system is supported by wind-driven upwelling, a bulk transport process that forces surface water offshore to be replaced with deep, nutrient-rich water 9,[20][21][22] .Upwelling favorable winds exhibit spatiotemporal variability, driving synoptic (shorter-term fluctuations lasting only for days or weeks) and seasonal cycles in upwelling, nutrient concentration, and oceanic phytoplankton density 9,23 .Additionally, global oceanographic and meteorological phenomena, such as the El Niño-Southern Oscillation and the Pacific Decadal Oscillation, alter productivity by introducing environmental variability on time scales of 4-7 and 20-30 years, respectively 24,25 .Finally, wind-driven upwelling is being modified by global climate change 18,26,27 , resulting in significant changes in coastal phytoplankton concentrations 28 ; further, pelagic-intertidal coupling can be expected to change with alterations in the winds, waves and stratification that control exchange between the coastal ocean and shoreline habitats.This complex array of interacting processes makes the variability of oceanic phytoplankton subsidizing intertidal communities extremely difficult to predict or model from first principles.Therefore, we propose the use of empirical oceanic phytoplankton time series data to represent the variability of the phytoplankton subsidizing intertidal food webs.Specifically, our present contribution uses a time series of oceanic phytoplankton (measured as chlorophyll-a) from the central coast of Chile to evaluate the effects of such variability on the dynamics of a rocky intertidal food web.
The bioenergetic model 29 expanded to large food webs by the Allometric Trophic Network (ATN) model 30 has been successful in capturing the dynamic processes of aquatic food webs, from predicting interaction strengths of rocky-intertidal food webs 31 to modeling the food web dynamics of lakes 32 .This model uses biomass (energy) as a currency and has been greatly influential toward understanding mechanisms behind community function and stability 33,34 and for predicting consequences of environmental changes and human exploitation on diversity 6,[35][36][37][38][39] .A main reason for its success is that the ATN model leverages several trophic and metabolic processes that scale with the body mass of aquatic organisms [40][41][42][43] , allowing researchers to use allometric scaling to parameterize the model for any system where body masses of the interacting species are known.This model, however, only accounts for the autonomous dynamics that emerge from allometrically scaled demographic and interaction rates.Previous work shows how the dynamics of the ATN model applied to the rocky intertidal system we study here are affected by including phytoplankton subsidy at a constant rate (fixed throughout the entire simulation), and subsequently increasing or decreasing this rate based on expectations from different climate change scenarios 6 .Here, we advance the field by modeling these subsidies more realistically by incorporating their high temporal variability, using time-series of chlorophyll-a to incorporate the empirical variability of oceanic phytoplankton into the subsidies received by the rocky intertidal food web.
Our specific objectives are to evaluate how the rocky intertidal food web model of the central coast of Chile responds to (1) the variability of pelagic phytoplankton abundance and (2) the rate of pelagic-intertidal mixing.These two factors combine to determine the flux of particulate food from the ocean to intertidal habitats.We found that while this pelagic subsidy is essential for fueling the intertidal food web, the high level of temporal variability (enhanced by rapid pelagic-nearshore exchange) increases biomass fluctuations putting species at increased risk of extirpation, and that variation in abundance due to phytoplankton subsidy is most pronounced in top predators relative to other types of species.

Results
Figure 1 outlines the study setup.Open-water chlorophyll-a measurements were obtained for the research site near Las Cruces, Chile (Fig. 1A).These empirical values were scaled to model units (biomass per area measured in g/m 2 , displayed as g/m 2 or kg/m 2 as denoted), and a spline curve was fit to the data (Fig. 1B) to create a continuous curve of offshore phytoplankton biomass.A network representation of the intertidal ecosystem is provided in Fig. 1C with nodes representing the intertidal species and links representing their associated trophic interactions.The height of a species' node corresponds to its trophic level.Variable pelagic phytoplankton (offshore phytoplankton, OP) drives this food web network by subsidizing phytoplankton in the intertidal habitat (foodweb phytoplankton, FP).The rate constant, k mixing , controls the rate at which pelagic phytoplankton is mixed into the food web network.k mixing is assigned values of 0.1, 1.0, and 10 h −1 as denoted.For ease of illustration, we will focus on results from 2003; however, trends are consistent across all years modeled.Results from other years are included in the supplementary materials and referenced where relevant.Summary information for the empirical dataset is provided in Supplementary Table S1.

Species' biomasses respond differently to phytoplankton variability depending on their food web connections
We began our investigation into the effects of pulses of high pelagic phytoplankton on intertidal population dynamics using an intermediate rate of pelagic-intertidal mixing, k mixing = 1.0 h −1 .
Pulses of high pelagic phytoplankton, such as the bloom event in late 2003 (Fig. 2A), are characteristic of seasonal upwelling in central Chile 23 .Intermediate to high peaks in offshore phytoplankton were observed in 10 of the 12 years modeled (Table 1, Supplementary Figs.S1, S2).We found that the abundances of all filter feeder species were strongly and positively influenced by increased phytoplankton subsidy (Supplementary Fig. S3).Our model assumes that intertidal particulate organic matter (comprised of phytoplankton and baseline detritus) is the only resource available to filter feeders.Consequently, the biomasses of filter feeder species are tightly coupled with the abundance of intertidal phytoplankton.Figure 2B shows the response of the filter-feeding barnacle, B. laevis, an important food source for omnivore sea snails, to the elevation in offshore phytoplankton that occurred in late 2003.
The total biomass of all producer species (including phytoplankton) remained relatively constant across the year (labeled Prod in Fig. 2C), but the total biomass of algal species was negatively influenced by increased phytoplankton subsidy (labeled Algae in Fig. 2C).These two results are explained by the community-level carrying capacity that is shared by all producer species of the model (see "Methods" section, Eq. 3).This carrying capacity limits the total biomass of producers such that when a producer species increases its biomass others will decrease theirs due to competition for space, light, or other resources.These results relating phytoplankton subsidy with total producer and non-planktonic producer biomass were consistent across all years modeled (Supplementary Fig. S1).
Responses among consumer species were more varied.For example, different sea snail species exhibited different responses to the Oct.-Nov.2003 elevation in phytoplankton based on their feeding preference as either omnivore or herbivore (Supplementary Fig. S4).On the one hand, omnivorous sea snails like L. orbignyi (LO in Fig. 2D) directly consume filter-feeding barnacles, placing them two nodes above food web phytoplankton.Consequently, our model shows biomass variations in omnivorous sea snails similar to those observed in phytoplankton and their filter feeder prey.On the other hand, herbivorous sea snails like T. atra (TA in Fig. 2D) consume algae but not barnacles nor any other consumers of phytoplankton.Consequently, increased phytoplankton affected herbivorous sea snails only indirectly by increasing the abundance of their predators and suppressing the growth of their resources.Through this example, we observe that offshore phytoplankton alters both top-down and bottom-up forces differentially even for related (within class) species, and that these differences can be made evident through inspection of their network connections.This trend-herbivorous sea snails decreasing their biomass in response to elevated offshore and foodweb phytoplankton and omnivorous sea snails increasing their biomass in response to elevated offshore and foodweb phytoplankton-is consistent across all twelve years modeled (Supplementary Fig. S1).

Increased pelagic-intertidal mixing increases biomass fluctuations putting species at increased risk of extirpation
We simulated the effects of varying the pelagic-intertidal mixing rate on the rocky intertidal ecosystem by assigning values of 0.1, 1.0, and 10 h −1 to k mixing (Fig. 3A,B).Although increasing values of k mixing increased the concentration of phytoplankton in the rocky intertidal (Supplementary Fig. S5), our simulations showed only a small relative increase in total primary productivity due to the community-level carrying capacity shared by producer species (Fig. 3C, Prod).However, total consumer biomass-and biomass of filter-feeding invertebrates in particular-increased with k mixing , highlighting the importance of phytoplankton in fueling the intertidal ecosystem across trophic levels (Fig. 3C, FF and Cons).
The pelagic-intertidal mixing rate exhibits profound influence on the persistence of certain species, as highlighted by the response of the sea snail Onchidella to peaks in offshore phytoplankton (Fig. 4, Supplementary Fig. S2).Using an intermediate mixing rate (k mixing = 1.0 h −1 ), our model predicts that the biomass of Onchidella would have reached a critical threshold in late October 2003 (Fig. 4B-solid line) caused by an increase in offshore phytoplankton (Fig. 4A) that produced a moderate reduction in Onchidella's resources (algae, Fig. 4C-solid line) and a large increase in its predator's biomass (Fig. 4D-solid line).The ATN model considers a species at risk of local extirpation when its calculated biomass drops below a critical threshold-in our case 10 -6 g/ m 26,34 -because such extremely low biomass levels do not recover without an external input of biomass of such species.Local extirpation of Onchidella was simulated at k mixing = 1.0 h −1 in six of the 12 years simulated (1999,  2000, 2002, 2003, 2006, and 2007).In all six of these years, offshore phytoplankton peaked at concentrations over 70,000 g/m 2 (Table 1, Supplementary Fig. S2).Decreasing the mixing rate to 0.1 h −1 results in slower exchange between intertidal habitat and offshore waters and rescues Onchidella from extirpation (Fig. 4B-dashed) because the biomasses of its resources (Fig. 4Cdashed) and predator (Fig. 4D-dashed) are significantly less affected by offshore phytoplankton.By lowering our simulated pelagic-intertidal mixing rate, Onchidella was rescued from local extirpation in all six of the years with high offshore phytoplankton (Table 1, Supplementary Fig. S2).Conversely, increasing the mixing rate to 10 h −1 results in almost immediate exchange with offshore waters and exacerbates the extirpation that was observed for intermediate mixing with both a stronger reduction in Onchidella's resources (Fig. 4C-dot-dashed) and a more dramatic increase in its predator biomass (Fig. 4D-dot-dashed).The biomass of Onchidella itself was much more variable under higher mixing rates with two additional incidents when Onchidella was at critical risk of extirpation (Fig. 4B-dot-dashed, triangles).Simulations performed at this elevated pelagic-intertidal mixing rate predicted a local extirpation of Onchidella in four additional years which showed only intermediate levels of offshore phytoplankton between 47,500 and 70,000 g/m 2 (2001, 2004, 2005, and 2008) (Table 1, Supplementary Fig. S2).
In our simulations, the survival of Onchidella was influenced by two factors: the peak abundance of offshore phytoplankton and the choice of pelagic-intertidal mixing rate.Our simulations indicate (1) that Onchidella can survive during periods with relatively modest peaks in offshore phytoplankton at all rates of pelagic-intertidal mixing, and (2) that Onchidella can survive scenarios with low pelagic-intertidal mixing rates regardless of the intensity of fluctuations in offshore phytoplankton abundance (Table 1, Supplementary Fig. S2).However, during years with intermediate peaks in offshore phytoplankton, Onchidella faced extirpation when subjected to high rates of pelagic-intertidal mixing (k mixing = 10 h −1 ).Further, during years with very high peaks in offshore phytoplankton abundance, Onchidella experienced extirpation at both high (k mixing = 10 h −1 ) and intermediate (k mixing = 1.0 h −1 ) levels of pelagic-intertidal mixing.This can be explained by the modelled effects of elevated foodweb phytoplankton has on summed algal biomass.Since both increased k mixing and elevated offshore phytoplankton act to increase foodweb phytoplankton concentrations, both factors will suppress the biomass of algal species.If the biomass of algal species drops too low, Onchidella is effectively resourceless and put at risk of extirpation (Fig. 3C).We observe that Onchidella's survival is modulated in a dose-dependent manner by peak abundance of offshore phytoplankton and that an appropriate combination of peak offshore phytoplankton abundance and pelagic-intertidal mixing rate is necessary for its survival.

Variation from phytoplankton subsidy is most pronounced in consumers and is exacerbated by higher pelagic-intertidal mixing rates
Within-year variations in species biomasses caused by fluctuations in offshore phytoplankton concentration differed across species guilds.To assess the variability across the network, each species was assigned to one of five guilds: algae, filter feeders, herbivores, omnivores, and carnivores.Variability within each guild was evaluated using a species-specific coefficient of variation that was normalized against the variability in offshore phytoplankton (Fig. 5).
For each rate of pelagic-intertidal mixing, the algal guild exhibited the lowest normalized coefficient of variation.The high biomass of algae and kelp causes the variation-to-mean ratio to be small relative to that of phytoplankton.Herbivore species also exhibited relatively low normalized coefficients of variation, explained by their tight trophic relationship with non-planktonic producer species.
Filter feeders, omnivores, and carnivores exhibited differential patterns of variability under different scenarios of pelagic-intertidal mixing.As the pelagic-intertidal mixing rate parameter was raised, variability in the biomasses of filter feeders increased.This can be attributed to two factors: (1) the strong trophic connection between filter feeders and the abundance of phytoplankton in the food web and (2) our previous result showing that higher values of pelagic-intertidal mixing rate enhance the availability of phytoplankton as a food source for filter feeders.Together, these two factors increase variability in filter feeder abundances, pushing variability closer to that of offshore phytoplankton.Variability arising from fluctuations in offshore phytoplankton abundance is more ambiguous in omnivore species.Omnivores exhibit extremely low relative variability under low pelagicintertidal mixing rates, with relative variability increasing dramatically with mixing rate.Variability in carnivore species' biomasses, however, shows elevated relative variability even in scenarios with low pelagic-intertidal www.nature.com/scientificreports/mixing rate.Further, carnivores show a clear trend of increased variability with increasing pelagic-intertidal mixing rate.Carnivores tend to have low total biomasses and incorporate several different sources of variability due to their position near the top of the food web, causing their biomass variability to be large.These observations are consistent with previous studies showing that predators amplify biomass fluctuations and exhibit high sensitivity to network perturbations 44 .

Discussion
The rich and diverse nearshore ecosystems of the rocky shore would not be possible without the subsidy provided by offshore phytoplankton 3,45 .However, the rate at which pelagic phytoplankton is delivered to intertidal habitats is highly variable, depending on both the concentration of offshore phytoplankton, which can vary by more than two orders of magnitude in upwelling systems, and the variability in pelagic-intertidal mixing rates, which are driven by a combination of waves, tides, and wind-driven currents.
Traditional ATN modeling using a stationary approach or a constant rate of subsidy does not capture the population dynamics that result from variable phytoplankton input to the network nor the food web level consequences of the different pelagic-intertidal mixing rates.Several recent publications have sought to incorporate time-dependent abiotic forcing into trophic network modeling using either wholly or partially simulated parameters including carrying capacity, intrinsic growth rate of producers, and metabolic rate of consumers among others [46][47][48] .Our work advances population dynamic modeling by developing a novel ATN methodology using an empirical time-series of chlorophyll-a to realistically simulate time-dependent subsidy (allochthonous input) in nearshore foodweb structure.This methodology allowed us to discover different responses to the variability in phytoplankton delivery that depends on the position of each species in the food web-and to relate those responses to oceanographic events reflected in the empirical time-series of chlorophyll-a.
A principal finding is that incorporating pulsed, time-dependent subsidy of phytoplankton into our population dynamics models influences transient food web dynamics and leads to important temporal fluctuation in species' biomasses.Moreover, our model suggests that the local rate of delivery of pelagic phytoplankton to the intertidal can amplify or alleviate the magnitude of these fluctuations, having important consequences on species' risk of extinction as illustrated by our results with Onchidella.High delivery rates essentially increase the magnitude of the subsidy pulse and thus amplify the magnitude of transient effects leading to greater temporal variability in populations.Our results reinforce those of others 44 that report carnivores, as a multispecies functional group, are particularly sensitive to perturbations, showing the overall greatest cumulative amplification of variation in biomass fluctuations due to elevated phytoplankton abundances.For intermediate species, this can lead to intensification and further amplification of top-down effects.For herbivore species not buffered by increased bottom-up effects, our simulations imply that the boost in phytoplankton pulse size can lead to sufficient intensification of top-down, negative effects that it pushes the population over a threshold of low biomass to a point that it cannot recover.These results are consistent with experimental observations showing increased risk of extirpation in meso-predator populations due to an imbalance in top-down and bottom-up forcing [49][50][51] .For these, variation in risk to extirpation depends on network connectivity and predator identity.In contrast, under similar pelagic phytoplankton abundance, simulations performed using lower mixing rates can buffer these populations from extinction by isolating them from offshore fluctuations and thus favor more balanced top-down/bottom-up effects that dampen biomass fluctuations.
Consistent with the work by Avila-Thieme et al. 6 using the same network model but with constant subsidy rates, we found that increased phytoplankton biomass decreases the biomass of non-planktonic producers via competition for the shared carrying capacity among producers.Additionally, both models found that increased phytoplankton biomass increased the biomass of filter feeders and that of their predators but decreased that of herbivores via the combined effects of increased predation (more abundant predators shared with filter feeders) and decreased resource availability (non-planktonic producers).Finally, both models found that omnivore species demonstrate greater resilience to elevated phytoplankton levels than herbivores by exploiting additional resources that are not affected by the decline in non-planktonic producers.This consistency between the two models suggests that these food web responses to increased phytoplankton subsidies are robust to the transient dynamics caused by the time-dependent pulsing of subsidies.Rather, these transient dynamics predict species extirpations, especially at higher trophic levels, that the model does not show using constant subsidy rates.We note these extirpations occur despite model baseline phytoplankton levels that support all species in steady state.Increased risk of extinction under transient dynamics is particularly important for understanding the effects that extreme oceanic events driven by climate change may cause to intertidal ecosystems.
No model completely represents all dynamics in a real system, and there are several points where additional empirical data may improve our simulations.The ATN model as used here only accounts for network structure, bioenergetic and trophic relationships between species, and empirical time-series data on chlorophyll-a concentration.Our choices of fluorescence scaling and pelagic-intertidal mixing rate rely on the previous model's subsidy value and may be improved using oceanographic simulations or empirical data.Without any available data recording the pelagic-intertidal mixing of nearshore phytoplankton, we estimate surf-zone flushing times from prior studies of circulation and dye concentrations 52,53 .We incorporate detritus through a baseline particulate organic matter concentration that supports filter feeders during extended periods of low offshore phytoplankton.While this dampens biomass fluctuations, they are sufficient to show significant results as demonstrated through our analysis of Onchidella.Additionally, there are several other biotic and abiotic factors that are not incorporated in our model-most importantly, sea surface temperatures, which modulate the biological rates underpinning ATN modeling, and nutrient concentration, which controls algal growth rates.
Further, to fully assess the vulnerability of food webs to time-varying subsidies, there is a need to account for the different time scales of specific trophic processes in relation to the time scales of environmental subsidies.For example, time-dependent changes in feeding and assimilation processes are not included in our model.Feeding and the capacity to capitalize on accessible prey can exhibit substantial variability owing to alternative saturation rates or shifts in behavior and foraging tactics aimed at striking a balance between risk and reward in the presence of multiple stressors.Additionally, differences in the timing of resource peaks may further contribute to the intricacies of this variability.Thus, future work to incorporate the relative timing of plankton subsidy to that of temporal variability of foraging could prove insightful to the sensitivity of propagating effects and whether timing-mismatches are enough to alter the relative balance of top-down and bottom-up effects.
Our study introduces a novel methodology incorporating empirical timeseries data into the ATN framework to study the effects of variable abundances of oceanic phytoplankton on species in a Chilean rocky intertidal ecosystem.Our results draw attention to the importance of transient population dynamics and the role they may play in predicting species extirpations during extreme events.Specifically, large fluctuations in offshore phytoplankton biomass put intertidal species at risk of local extirpation even though levels of intertidal phytoplankton were not allowed to drop below a level that supported the persistence of all species in steady-state simulations.Our results suggest that for a species to persist in the face of large-scale biomass variability introduced by bottomup processes, that a concomitant top-down, network-scale process is required to dampen biomass oscillations allowing for species persistence.

Chlorophyll-a dataset and preparation
Nearshore chlorophyll-a density ( δ chloro g/m 3 ) was measured daily from a permanent sampling quadrat at the Estación Costera de Investigaciones Marinas (ECIM) in Las Cruces, Chile between Jan 1999 and Dec 2010 28,54 .We applied a linear approximation with slope β 1 and intercept β 0 to transform chlorophyll-a density to phytoplankton population density ( D phyto individuals/m 3 ).We used a frequency-weighted average body size ( m phyto g/individual) of several phytoplankton species present at the Las Cruces sampling site 6 to convert from phytoplankton density to volumetric phytoplankton density ( δ phyto g/m 3 ) according to We normalized the total amount of phytoplankton subsidies to the time-variable model over the twelve-year study period to the amount of subsidies that would be delivered to the Avila-Thieme 2021 model over the same time interval 6 .Avila-Thieme 2021 uses a daily phytoplankton subsidy of 7355 g m 2 d .The Avila-Thieme 2021 model delivers a total phytoplankton subsidy, s total , to a quadrat with area A m 2 of The time-variable model provides phytoplankton subsidy to the volume under a quadrat with volume V m 3 and average depth z ave .Assuming the foodweb compartment is rapidly depleted of chlorophyll to baseline, we calculate a total phytoplankton subsidy to the time-variable model of where α is a new scaling factor representing three distinct physical properties: the mass of individual phytoplank- ton, the scaling from chlorophyll density to phytoplankton density, and the average depth of the water column under the quadrat.α has units of m 1 according to Setting the total phytoplankton subsidies equal between the two models, we obtain a scaling factor of α = 1.567m.We scaled our empirical chlorophyll density time-series (g/m 3 ) by alpha (m) to obtain an area- based approximation of phytoplankton biomass ( B phyto g/m 2 ).A continuous spline was fit to the phytoplankton biomass time-series and was numerically differentiated to obtain the rate law governing the biomass of offshore phytoplankton, B op (see Network Structure Changes under Non-Autonomous ATN Extension below).

Allometric trophic network
The Las Cruces food web contains 106 species with 1,362 trophic links (Fig. 1A, ref. 55 ).This foodweb includes all species found to co-occur during community structure surveys carried out at multiple, wave-exposed rocky intertidal sites along a ~ 700 km stretch of the central Chilean coast 16,56,57 .All 106 species are known to coexist at the Las Cruces site 6,57 .
. www.nature.com/scientificreports/Our ATN modeling utilizes two sets of differential equations to describe the biomass of each producer (Eq. 1) and consumer (Eq.2) species (recorded in units of grams biomass per square meter) according to the principles of mass balance as follows: The biomass of producer species changes according to the balance of autotrophic growth with predation (Eq.1).Autotrophic growth is controlled by the intrinsic growth rate, r i , and current biomass, B i , and is scaled by a logistic growth factor, G i (Eq.3).The logistic growth factor takes into consideration the shared, communitylevel, carrying capacity for producer species, K, to slow producer growth as the summed producer biomass approaches or exceeds carrying capacity.Biomass loss to predator j from prey i increases with the predator's mass-specific metabolic rate, x j , and prey-specific attack rate, y ji , and decreases with the assimilation efficiency of prey j into predator i, e ji .
The biomass of consumer species changes according to the balance of consumer growth rate, explicitly modeled metabolic loss and mortality, and loss to predation.Consumer growth rate is modeled as assimilation of prey biomass; consequently, consumer growth rates are a function of their total consumption rate, x i B i jǫprey y ij F ij (B) , scaled by an assimilation coefficient, f a .Metabolic loss including mortality is directly pro- portional to current biomass and metabolic rate and is scaled by a mortality coefficient, f m .Finally, the loss to predation term is identical to that of producers.The Holling's functional response, F ij , is used to determine the consumption rate of each resource for each consumer (Eq.4).This function considers the relative preference of consumer i for resource j, ω ij , as well as the Holling's coefficient, q, to determine curve shape.We invoke an intermediate Holling's Type II and Type III response by using q = 1.2 58 .
Supplementary Table S3 shows all model parameters with their initial values and equilibrium values.The Las Cruces ATN model was parameterized as in ref. 6. Intrinsic growth rate of producers, metabolic rate of consumers and maximum consumption rate of consumers were parameterized using measures of each species' body size 6,55 .Initial biomass values were estimated using species density (mobile species and cnidaria) and surface cover (sessile species) from six years of empirical records maintained at the research site 6 .Community-level carrying capacity was estimated as the sum of biomasses for the fastest growing species in each of six producer functional groups (microalgae, ephemeral, corticate, crustose, coralline, and kelp) weighted by the number of species in each group 6,59 .We investigated the effects of varying the community-level carrying capacity on community composition (Supplementary Fig. S6).Half-saturation density parameters were estimated from past literature using the Lake Constance ATN model as a primary reference 6,32,60,61 .

Non-autonomous ATN extension
Parameter Modifications: We forced our ATN model to use a timestep equal to one hour.Consequently, we converted intrinsic growth rates of all producers, r i , metabolic rates of all consumers, x i , and maximum consumption rates of consumers, y i , to an hourly rate using the daily rates reported in ref. 6 .All other parameters, including initial conditions, are identical to those used in ref. 6 .
Network Structure Modifications: To encapsulate the phytoplankton dynamics of this system, we split the existing phytoplankton node into three nodes: offshore phytoplankton, foodweb phytoplankton, and "baseline phytoplankton" representing detritus.
1.The biomass for offshore phytoplankton, B op , is equal to the scaled, empirical dataset described in Chloro- phyll-a Dataset and Preparation.This node connects to the food web only by water-borne exchange with the foodweb phytoplankton node.This introduces a new rate constant, k mixing , which controls the rate at which plankton moves from offshore into foodweb phytoplankton.k mixing is assigned values of 0.1, 1.0, and 10 h −1 as denoted.2. We separated the existing phytoplankton node into two nodes, baseline phytoplankton and food web phytoplankton.These two nodes share identical connections within the food web. (1) Loss to predation j x j y ji B j F ji (B) e ji .
( www.nature.com/scientificreports/ The biomass of "baseline phytoplankton" (representing detritus contribution to particulate organic matter available to filter feeders) is held constant at a value of 3,750 g/m 2 .This node supplies the minimal amount of phytoplankton required to prevent network collapse under steady state conditions (Supplementary Figs.S7, S8).This node represents organic particulate matter available for consumption by filter feeders such as detritus that is otherwise not modeled in our ATN.
The foodweb phytoplankton node provides the primary supply of phytoplankton to consumers in the food web (Supplementary Fig. S5).This node encapsulates the local dynamics of intertidal phytoplankton.Further, this node is subsidized by biomass from offshore phytoplankton via water-borne exchange (Eq.5).This subsidy modifies the ATN foodweb phytoplankton according to the following relation.

Simulations
Each year of empirical data, from Jan 1 to Dec 31, was analyzed independently.To better compare simulations performed under different empirical datasets, we sought to reduce the transient dynamics introduced by having initial conditions that were significantly different from runtime conditions.We evaluated two separate strategies for reducing the impacts of initial conditions and found little change between the two.Strategy 1 involved running simulations under an "equilibration phase" with the offshore phytoplankton node deactivated for ten years of simulation time prior to activating the empirical data during the "experimental phase" (Supplementary Fig. S9).Strategy 2 replaces the "equilibration phase" with ten appended copies of the empirical dataset.This allows population biomasses to approach runtime conditions more gently than strategy 1. Population distributions were consistent between each strategy, and the results presented here were obtained using strategy 1.

Coefficient of variation and violin plots
We use the coefficient of variation, the ratio of standard deviation to mean, as the preferred metric of variability for the timeseries data presented here.We calculated coefficient of variation for each species during the time interval when its biomass is above the extinction threshold of 10 −6 g/m 2 (i.e., trailing zeroes were truncated where applicable), and normalized these values to the coefficient of variation that was calculated for the biomass of offshore phytoplankton.
Violin plots were used to display distributional data on coefficients of variation.These plots are a modification of box plots in that they show the median value, each individual value, interquartile range, and shade the kernel density estimate to give an indication of the data distribution.Wider sections of violin plots represent a higher probability that members in the population will take on a given value.Violin plots were generated in Matlab using code from ref. 62

Figure 1 .
Figure 1.(A) Empirical timeseries dataset of offshore chlorophyll-a measurements in 2003.(B) Continuous spline fit to scaled empirical dataset.(C) Graphical representation of the Las Cruces food web network.Nodes are vertically arranged by using prey-averaged trophic level.Important nodes are labeled as follows: OP offshore phytoplankton, FP food web phytoplankton, BL the barnacle, B. laevis, Onc the sea snail, Onchidella sp., TA the sea snail, T. atra, LO the sea snail, L. orbignyi.

Figure 2 .
Figure 2. Biomass curves simulated using empirical data from 2003 for (A) offshore phytoplankton (OP, solid line) and food web phytoplankton (FP, dashed line), (B) the barnacle, B. laevis (BL), (C) summed biomass of all producer species (Prod, dashed line) and algal species only (Algae, solid line), and (D) the omnivorous sea snail, L. orbignyi (LO), and the herbivorous sea snail, T. atra (TA).Biomasses are shown in g/m 2 .These results use an intermediate pelagic-intertidal mixing rate, k mixing = 1.0 h −1 .

Figure 3 .
Figure 3. (A) Measured biomass of offshore phytoplankton from Sept 2003 to Jan 2004.(B) Simulated biomass of food web phytoplankton under k mixing = 0.1 h −1 (dashed line), 1.0 h −1 (solid line), and 10 h −1 (dot-dashed line).(C) Range of simulated biomasses from Sept 2003 to Jan 2004 for four trophic categories and three values of k mixing .This range spans the minimum simulated biomass to the maximum simulated biomass.Prod summed biomass of 46 producer species, FP food web phytoplankton biomass, Algae summed biomass of 45 algal species, FF summed biomass of 15 filter feeder species, Cons summed biomass of 60 consumer species.All biomasses are reported in kg/m 2 .

Figure 5 .
Figure 5. Violin plots showing the coefficients of variation for each species normalized to that of offshore phytoplankton.Results shown for 2003 and each value of k mixing investigated.Species were grouped into one of five guilds: algae (Algae), filter feeders (FF), herbivores (Herb), omnivores (Omni), and carnivores (Carn).
fp G fp (B) − � j x j y j,fp B j F j,fp (B) e j,fp+ k mixing � B op − B fp � , B op ≥ B fp r fp B fp G fp (B) − � jx j y j,fp B j F j,fp (B) e j,fp , B op < B fp .

Table 1 .
Results list peak offshore phytoplankton abundance and simulated extirpation status of Onchidella under each pelagic-intertidal mixing scenario investigated.Cells with an 'e' signify simulations in which Onchidella experienced extirpation, and empty cells signify simulations where Onchidella survived.